function M = mass(dofs,mesh);
%        M = mass(dofs,V,T,posT,d,areas);
% this is for single degree version;
% % get the areas of all triangles.
% x21 = V(T(1:posT,2),1) - V(T(1:posT,1),1);
% x31 = V(T(1:posT,3),1) - V(T(1:posT,1),1);
% y21 = V(T(1:posT,2),2) - V(T(1:posT,1),2);
% y31 = V(T(1:posT,3),2) - V(T(1:posT,1),2);
% areas = (x21.*y31 - x31.*y21)/2;
V=mesh.V;
T=mesh.T;
posT=size(T,1);
d=mesh.d;
areas=mesh.areas;
mat = build(d);  [i,j,s] = find(mat);  L = length(i);
m = (d+1)*(d+2)/2;  m = m*m*posT; 
S = zeros(m,1);
Indx1 = zeros(m,1); Indx2 = zeros(m,1); 
pos = 1; % it may be accerate by cancel the loop.
for tri = 1:posT    
    Indx1(pos:(pos + L-1)) = dofs(tri,i);
    Indx2(pos:(pos + L-1)) = dofs(tri,j);
    S(pos:(pos + L-1)) = areas(tri)*s;
    pos = pos + L;    
end
dim = max(max(dofs));
M = sparse(Indx1(1:(pos-1)),Indx2(1:(pos-1)),S(1:(pos-1)),dim,dim);